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ABSTRACT 

An approximate technique is presented for the prediction of 
the large-scale turbulent structure sound source in a super- 
sonic jet. A linearized Euler equations code is used to solve 
for the flow disturbances within and near a jet with a given 
mean flow. Assuming a normal mode composition for the 
wave-like disturbances, the linear radial profiles are used in 
an integration of the Navier-Stokes equations. This results 
in a set of ordinary differential equations representing the 
weakly nonlinear self-interactions of the modes along with 
their interaction with the mean flow. Solutions are then 
used to correct the amplitude of the disturbances that rep- 
resent the source of large-scale turbulent structure sound in 
the jet. 

INTRODUCTION 

While considerable progress has been made over the last 
decade toward the prediction of jet noise, an efficient tool 
for the prediction of jet noise in practical geometries is 
still lacking. The problem can be split into the predic- 
tion of the source field and the prediction of the acoustic 
field associated with it. As for the latter, there are now a 
number of schemes that can successfully predict the acous- 
tic field associated with a well-prescribed sound source. 
In the acoustic analogy approach, time-average correla- 
tions can be used to construct the acoustic source. 1 This 
source is used in a volume integration to obtain the acous- 
tic field. An implementation of this approach is described 
by Khavaran 2 where the statistical properties of the noise 
sources are described by two-point, space-time correlation 
models derived from RANS solutions for the flow field us- 
ing k-e turbulence models. An alternative approach was 
recently developed by Tam & Auriault 3 where models for 
the two-point correlations are given in a fixed reference 
frame, rather than the moving reference frame commonly 
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used in previous approaches. In both of these approaches, 
it is the fine-scale turbulence that is modeled as the source 
of noise. 

Strong evidence suggests that jet noise, particularly in 
the supersonic regime, has contributions from large-scale 
turbulent or wave-like structures in the initial region in the 
jet and can be the dominant noise source. These struc- 
tures cannot be captured by classical turbulence modeling. 
Direct Numerical Simulations (DNS) and Large Eddy Sim- 
ulations (LES) can successfully capture these structures, 
but they are computationally intensive. 

Large-Eddy Simulations (LES) is the most promising 
approach for prediction of the noise source. The resolu- 
tion in LES involving the acoustic field is usually selected 
such that all the acoustically relevant scales are resolved. 
In principal, LES can be extended to resolve all the scales 
as in DNS by making the resolution equivalent to that of 
DNS. At the other end, the computations can be made less 
CPU intensive by selecting a coarse mesh, and resolving 
only the very-large scales (VLES). Although VLES may 
be less CPU intensive compared to LES, still there is an 
obvious need to develop an approximate, but fast approach 
for prediction of jet noise. 

The work presented herein focuses on predicting the 
noise source associated with the large-scale, low frequency 
noise. This work, along with a parallel work for prediction 
of the small-scale noise, should provide a practical tool for 
prediction of jet noise associated with both the small and 
large scales. The basic idea herein is to first neglect non- 
linearity and use a Linearized Euler Equations (LEE) code 
to predict the time-dependent source in a given mean flow. 
Nonlinear effects are then accounted for via a nonlinear in- 
tegral technique similar to that of Dahl & Mankbadi. 4 

The starting point of this analysis is that we consider 
a turbulent round jet at a sufficiently high speed so that 
the compressibility is significant. The development of this 
jet in the unexcited case is assumed to be given by some 
other means (e.g. analytically or via Reynolds averaged 
numerical simulations). This jet is then excited by a single- 
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frequency instability wave. The nonlinear development of 
this wave will be presented herein based on the integral 
energy approach. Along with the wave development, the 
mean flow-spreading rate is also modified. Since the fo- 
cus here is the supersonic jet, we must consider the helical 
modes as they are more amplified than the axisymmetric 
ones. Once the computations are completed, a description 
of all the variables of the large-structure disturbance field in 
and near the jet is given for a particular mode number and 
frequency. It can be used as the acoustic source to compute 
the noise radiation field outside of the jet. This is described 
in a separate paper. 5 

In the next section, we describe the LEE method used to 
compute the large-scale turbulent disturbances. The advan- 
tage of this approach is that it allows the base mean flow to 
be non-parallel. This has the consequence that the numer- 
ical singularities, present in the previous locally-parallel 
stability analysis, 4 are eliminated. The development of the 
integral energy equations is described in the following two 
sections. Finally, results are shown for the nonlinear calcu- 
lations showing the effects of dissipation, initial large-scale 
structure amplitude, and jet Mach number on the spatial 
evolution of the amplitude and phase of the structure. 

LINEARIZED EULER EQUATIONS 

Consider a high-Reynolds number turbulent jet issuing 
from a nozzle of diameter D in a still air. The jet is 
shock-free, but the Mach number is high enough for com- 
pressibility effects to be significant. The density and the 
component velocities are normalized by the jet exit density 
and axial velocity at the centerline, pj and Uj , respectively. 
The pressure is normalized by PjVj, time by D/Uj , and 
spatial coordinates by D. Each flow parameter is split 
into a time-averaged part Ui(x, r, (f>) and a disturbance part, 
rt^( x, r, (f>, t). Thus, the velocity can be written : 

Ui = Ui(x,r,<j>) + Ui(x,r,<j>,t) (1) 

where i = 1,2,3. In the cylindrical coordinates, 1 refers 
to the axial direction x with axial velocity u, 2 refers to the 
radial direction r with radial velocity v, and 3 refers to the 
azimuthal direction <J> with azimuthal velocity w. An over 
bar, ( ), denotes a time-averaged quantity. The pressure 

and the density are similarly split: 

P = P{x,r,<f>) + p'(x,r,</>,t) (2) 

p = p(x,r,(j))+p'(x,r,<f>,t) (3) 

For the product of the density with the velocities, we get 

put = (p + p'){Ui + u'i) = pUi + pu'i + p' Ui + p'u'i 
Time-averaging yields 

pUi = pUi + p'u\ (4) 


and we define 

put = pui - pUi = pu' i + p 1 Ui + p'u'i - pfu'j. (5) 
Thus, we get 

pui =pu~i+pui. (6) 

Starting from the full Navier-Stokes equations in conser- 
vative form, we neglect viscosity, linearize about a given 
mean flow and separate the disturbance variables into az- 
imuthal modes by assuming the presence of an exp(in^) 
factor in each variable. The resulting Linearized Euler 
equations may be written in cylindrical coordinates as fol- 
lows: 


(rQ) t + ( rF) x + ( rG) r + inH = S 


(7) 


where n is the azimuthal mode number and the vector quan- 
tities are: 


F = 


Q = [p r pv pw pu e'] 
pu 

—p'V U + pvU + puV 
—p'W U + pwU + puW 
-p'UU + 2puU + p' 

U(p' + e') + (j5u - p'U) ( P±M 


G = 


pv 

-p'VV + ZpOV + p' 
—p'V W + pwV + pvW 
-p'VU + puV + pvU 

V(p' + e!) + (pw - p'V) ( P + E 


H = 


pw 

— p'V W + pi’ W + pwV 
-p'WW +2ptiW + P ' 
-p'UW + puW + pwU 


W(p' + e!) + ( pw - p‘ 


’w) 


S = 


0 


-p'WW + 2ptiW + p' 
p'VW -pwV -pvW 
0 
0 


where the mean and the disturbance pressures are given by 

P=( 7-1) 


E - ip([/ 2 + V 2 + W 2 ) 


P' = { 7 - 1 ) 


e' — (pu U + pvV + pw W) 


+ \p\U 2 + y 2 + W 2 ) 
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E is the mean total energy and e' is the energy disturbance 
variable. In developing the Linearized Euler equations, 
only the linear terms of equation (5) for pu,, are retained. 

In this formulation, each additional mode calculated 
adds another set of equations to be solved. The advantages 
to this formulation are lowered storage and computation, 
improved centerline behavior, and improved boundary con- 
dition specification. 

Numerical algorithm 

The code uses an explicit time-marching method coupled 
with central differences in space. To eliminate spurious 
high-wavenumber oscillations, artificial dissipation is em- 
ployed. For computational efficiency and maximum com- 
piler optimization, the code is written in FORTRAN 77. 

The time marching method is a low-storage, fourth-order 
extension of an optimized LDDRK 5-6 method . 6-7 For ac- 
curate spatial differencing, the 7-point DRP method of Tam 
& Webb is used . 8 For artificial viscosity, an explicit 10th 
order dissipation is used . 9 

Computational grid 

The grid used for the cases computed for this paper 
had 171 (radial) x 251 (axial) grid points. In the radial 
direction, the grid started with the minimum spacing of 
A r/D = 0.02 at the centerline and smoothly stretching 
to 0.15 at r/D = 16. In the axial direction, the minimum 
spacing was Ax/D = 0.02 at the x = 0 boundary and 
smoothly stretching to 0.15 at the x/D = 35 boundary. 
The maximum spacing corresponds to 10 points per wave- 
length, which is well within the accuracy range of this code. 
All cases used the same computational grid. A Courant 
number of 1.4 was used for these computations. 

Boundary conditions 

Special attention is given to the boundary treatment in or- 
der to avoid non-physical oscillations, which can render the 
computed unsteady solution unacceptable. Several bound- 
ary treatments were considered. The boundary treatments 
follow that of Mankbadi et al . 10 At the inflow, Thompson, 
non-reflecting boundary treatment is implemented. The 
conventional acoustic radiation condition is applied at the 
radiation boundaries, which are defined at inflow, x = 0 , 
where r/D > 2, and at the outflow, x = x max , for all ra- 
dial points where the Mach number is less than 0.01. For 
r = r max , the radiation condition is applied to all appoints. 
At the outflow, where the Mach number is greater than 
0.01, the Tam & Webb 8 asymptotic boundary treatment is 
applied. In this code, the centerline boundary is represented 
with a point at the centerline, and a ghost point reflected 
across the centerline in the radial direction. Without az- 
imuthal mode decomposition, the centerline treatment for a 
three- dimensional problem is not straightforward, and was 
addressed by Shih et al." However, using the azimuthal 
mode decomposition method, the centerline boundary con- 


dition becomes straightforward: 

r [Gr G 2 G 3 G 4 G 5 ]^ =0 = 0 

and 

r [Gi G 2 G 3 G 4 G 5 ;';_ a ,. 

= r[G 1 — G 2 — G 3 G 4 G b ] ^ =Ar e imr 

where G'i to G 5 represent the five terms in the G vector in 
equation (7). 

Inflow disturbance 

At the inflow boundary, x/D = 0, a disturbance based 
on the linear stability eigenmodes for the inflow jet profile 
is introduced: 

Q = L Q (r)e iu;t 

where Lq (r) represents the radial eigenvector of the linear 
stability wave for the set of disturbance variables. To in- 
troduce the input disturbance into the flow field, the time 
derivatives of the disturbance are added to the computed 
flow variables at each time step: 

Qt, final = Qt.BC + iuL Q (r)e iuJt 

This boundary condition inputs the loc ally-parallel eigen- 
function solution that is not exactly equivalent to the Euler 
solution. This difference results in a transient wave near 
the inflow boundary causing inaccuracies in the numerical 
results. 

WEAKLY NONLINEAR EFFECTS 

The source predicted by LEE neglects nonlinear effects. 
However, nonlinear effects seem to be important and need 
to be taken into account. Dahl & Mankbadi 4 developed a 
nonlinear technique for the prediction of the noise source 
associated with instability waves in a compressible jet. 
In their work, the partial differential equations were in- 
tegrated across the transverse direction assuming that the 
disturbance profiles were given by the eigenfunctions of a 
locally-parallel, viscous, linear stability calculation. Cer- 
tain mathematical singularities that occurred within the sta- 
bility equations of this approach limited its usefulness. In 
addition, the integral energy equations of the previous pa- 
per contain some inconsistencies that are resolved in the 
presentation below. In the present work, the disturbance 
variable profiles will be taken from the linear Euler solu- 
tion. 

Shape assumptions 

In the integral energy method, the system of equations is 
integrated in the radial direction using shape assumptions 
for the disturbance variables. The disturbances, considered 
to be the large-scale turbulent structures, are assumed to 
have the form of a traveling wave. Following the work of 
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Lee & Liu, 12 the waves are assumed separable into an un- 
known amplitude function and a radial shape function 

/ u'i(x,r, <j>,t)\ / Ui(r)\ 

p’(x,r,(j),t) = A(x)e l ^ x ^ p(r ) exp('I') (8) 

\p'{x,r,(/),t)J \p(r)J 


where the Reynolds number. Re = pjUjD/p, , is initially 
assumed to be constant. The subscripts denote differentia- 
tion, the Laplacian is 

<9 2 ld_ d_ 1 d 2 

dx 2 r dr dr r 2 90 2 


'T = — icot + iruj) 

with an axial phase function ip(x) that also needs to be 
determined. In equation (8), (") denotes the radial shape 
function of the transverse coordinate r at a given location 
along the jet. Ui(r), p(r), and p(r) are radial shapes com- 
puted by the LEE code at a given axial location and at a 
given n and u>. Here, n is the azimuthal wave number indi- 
cating the rotation around the jet centerline. A(x) is the real 
amplitude function of x and is to be determined by a non- 
linear analysis. For the LEE results, the variables contain 
both axial and radial information. Only radial information 
is desired for the integral technique. The radial shape func- 
tions are estimated at each axial location by normalizing in 
the radial direction 

(|ti| 2 + |u| 2 + |tu| 2 ) rdr = 1 (9) 

at each axial location. 

Equations of motion 

The formulation begins with the following nondimen- 
sionalized continuity and momentum equations in cylindri- 
cal coordinates: 

Pt + ( pu) x + ^ [r(pv)\ r + ^(pw)<p = 0 ( 10 ) 



(pu) t + (pu 2 +p) x 

+ - [r (P uv )] r + \ (P UW )<j> 

= itc Ali + * V ^ x 


(ID 


(pv)t + ( puv) x 

+ - [r (p + pv 2 )] r + 1 (pvw), - pw 


r 

■p + ± 

r Re 


An j ( v + 2u V) 


( 12 ) 




(pu’)t + (puw) x 


(13) 


1 


+ - [r (pvw)] r + -(p + pw 2 ), + 


1 

Re 


Ain j ( w ~ 2v <i>) 


pwv 
r 
1 1 
3Re r 


(V-v) 


and 


_ du 1 d 
V • v = — + -~^~rv - 
dx r dr 


1 dw 
r dcj) 


Making substitutions from relations (1) to (6) into (10) 
to (13), we get: 


(P + P') t + ( P u + P u ) x 


(14) 


+ - \r {pv + pv)] r + -(pw + pw ) <p = 0 


{pu + pu) t + [P + p' + (U + u') {pu + pu)] x 
+ - [r (V + «') (pu + pu)] r 
+ ~[( W + W ') (P a + pu)] (t> (15) 

= st A ( E7+ "') + i (v ' F+v '»' ) ' 


{pv + pv) t +[{u + u') {pv + pv)] X 

+ ~[ r ( P +P'+( V + V ') (P® + pv))\ r 

+ - \_{W + w') (pv + pv)] ^ — - (W + w') (pw + pw) 


= l(F + P') 


(16) 


1 

Re 


( V + v')-±(v + v' + 2(W + w '),) 


1 


+ — (V-F + V..'), 


(pw + pw) t + [( U + u) (pw + pw )] , 


-I — [r (V + v') (pw + pw)] r + - (V + v') (pw + pw) 


1 


H — \P + p' + (W + w') (pw + pw) 


(17) 


1 


1 


Re [A (W + w') - ^ [W + w 1 - 2 (V + v ’) , 

_LI (v .f + v-„',* 


These equations will be used to derive a set of equations 
governing the mean flow and a set of equations governing 
the large-scale structures. 
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Mean flow 

The mean flow equations are obtained by time-averaging 
equations (14) to (17) containing the two-component de- 
composition. The continuity equation for the mean flow 
is 

(W) x + ^[r C pv)] r + \ (P™)# = 0 (18) 

and the three momentum equations for the mean flow are: 


r (pu V + v'pu'j 


( puU + u'pu + p] H — 

V ) x r l 

+ i {flaw + w'piL S j = ^A(7 + — !^(V • V) x 


(19) 


Re 


3Re 


(pv U + u'pv'j 4 — r (pv V + v'pv + P^j 
4 — {flu W + w'pv^j (flw W + w'pv^j 


= -P+ — 
r Re 




( 20 ) 




(pw U + u'pw'j 4 — r (pw V + v'pw ^ 


( 21 ) 


4 — (pw W + w' pw + p\ + - (pw V + v'pw 
r \ J <b r \ 


1 

Re 


AW - 4t (W - 21.,.) 




Large-scale structure 


The equations governing the large-scale structures are 
obtained by subtracting the mean flow equations, (18) to 
(21), from equations (14) to (17). The continuity and mo- 
mentum equations for the large-scale structures are: 

p't + pu x + - ( rpv) r + - pw $ = 0 ( 22 ) 


flu t 4- (p' + flu v! + pu U + pu u! — flu u! 

+ - r (flu v' + pu V + puv' — puv , 'j 

H — ( ~puw ’ + puW + puw' — lmw'\ 
r V / <f> 

= s-Aw' + (V • v') x (23) 

Re 3Re v ,x 


pv t + (pv v! + pvU + pv v! — pv u'J 

r (p 1 +flvv' + pvV + pvv' — pvv’ 


r 


+ - 
r 


- (pv w' + pvW + pv w' — pu w ,S j 


Al/ -Z2( V ' + H) 


- ( pw w' + pw W + pw w' — pw w' 
r \ 

1 , 1 

: ~P + w~ 
r Re 


pw t + (pw u' + pw U + pw u' — pw u' 

4 — r (flw v' 4- pw V 4- pw v' — pw t>'^ 

1 ( , , — ~ , — 

4 — [p + pw w 4- pw W 4- pw w — pw w' 
r V 

1 / . _ — _ . 

-I — pw v + pw V + pw v — pw v' 
r \ 


( 24 ) 


1 


Aw ' ~ 02 ( W '~ 2 V 4>) 


1 

Re 




(25) 


Kinetic energy equations 

Mean flow kinetic energy equation 

The mean flow kinetic energy equation is obtained by 
first multiplying the x-momentum equation (19) by U, the 
r- momentum equation (20) by V, and the ^-momentum 
equation (21) by W. Then, the resulting equations are 

added together. After much manipulation, using the mean 

— 2 — 2 

flow continuity equation, and defining K = (U + V + 

■£ 2 . 

W )/2, the combined equation is rearranged to obtain the 
kinetic energy equation for the mean flow 

(flu K + puu' U 4 -puv'V + puw'W^j 


id_ 

r dr 


r (pv I\ 4- pv u' U 4 - pw' V 4- pv w' w\ 

1 d ( — — 2 = — — \ 

4 — — pw K + pw u' U + pwv' V 4- pw w' W 
r dtp V / 

W— 

4 -UP x + VP r + —P < t, 


— U x puu' — V x puv' — W x puw' — U r pvu' 

1 ^ 2 = 2 -^ 1 ^ 


— V rpv v' — W r pvw' U ppwu' V ppwv' 

r r 


1, 


R 


W ppw w' V pw w' 4 — W pw v' 

r r r 


1 

Re 


AK - (U ix ) 2 - (U ir f - 1 


( 26 ) 


4- 


3Re 


W . 


U{V ■ V) x + R(V • V) r + — (V • V)* 
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where ( U i X ) 2 = ( U x ) 2 + ( V x ) 2 + (VF^) 2 and similarly for 
( U ir ) 2 and 

Large-scale kinetic energy equation 

The large-scale kinetic energy equation is often obtained 
by multiplying each disturbance momentum equation by 
its corresponding disturbance velocity and then adding the 
three equations. This process was used when shape as- 
sumptions were based on using locally-parallel, stability 
calculation results. With LEE, we obtain disturbance vari- 
ables that are complex. In anticipation of using these 
complex values, we develop the large-scale kinetic energy 
equation by multiplying each disturbance momentum equa- 
tion by the corresponding conjugate disturbance velocity 
and then adding the conjugate, that is, for the x-momentum 
equation: 

«'*(x-mom.) +a'(a:-mom.)* (27) 

Similar equations are derived using the r-momentum equa- 
tion and the ^-momentum equation. After adding these 
equations together, manipulating extensively, and applying 
a time-average, the resulting kinetic energy equation for the 
large-scale structures is 


d 1 d \ d 

— (ZpuQ) + ~—(2rpvQ) + ~ — (2pwQ) 
Ox r or r dtp 


+ u'*p' x + v'*p' r + 


w'% 


;=^ — — — u'* pw — 

+ u'*pu U x + u'*pv U r -) U ft 

r 


— — — 1 =^ — v pw — 

+ v'*puV x + v'*pvV r H V a, 

r 


— ;=zr — w'* pw — 

+ w'*puW x +w'*pvW r H W $ 


w'*pw— v'*pw==- 
V H W + cc+ ht 


1 

Re 


A2Q - u'* x u' ix - u'*u' ir - 


U 'i4> U 'i<\) 


v'* w'* 

- 72 K + 2w '<t>) - - 2v 'd>) + cc 


(28) 


into two terms; one for the derivative of the amplitude func- 
tion A(x) and one for the derivative of the phase function 
ijj(i r). Thus, an additional equation is required to describe 
the axial evolution of the phase. Following Lee & Liu, 12 
this is obtained by first subtracting the momentum equa- 
tions, rather than adding as in equation (27), 

«'*(x-mom.) - *t'(x-mom.)* (29) 

This equation and the additional two equations for the r- 
and </>-momentum equations are added together and exten- 
sively manipulated. The resulting equation is more com- 
plicated than equation (28) and will not be given. Only the 
final integrated equation will be given below. 


INTEGRAL FORM OF THE ENERGY EQUATION 

For the round jet, the mean quantities are assumed to 
be axisymmetric. Thus, W = 0 and <9( )/d<f> = 0. The 
energy equations (26) and (28) for the mean flow and the 
large-scale structures, respectively, are simplified as these 
terms are removed. Then, the usual boundary-layer-type 
approximations are applied to the mean quantities. The 
mean radial velocity is much less than the mean axial ve- 
locity, V <C U, and the axial gradients are much less than 
the radial gradients, d( )/dx <C d( ) /dr. Using these ap- 
proximations and considering the number of disturbance 
variables in each term, the terms are ordered in size and 
those terms up to second order are retained. The mean flow 
pressure is then assumed to be constant across the jet. Fi- 
nally, multiplying by r and integrating over r, the integral 
form of the mean energy equation reduces to the following 
simple form. 


d f°° ( 1 — 2 . , 

-puU + puu' U I rdr 


dx 


to 


= / puu 1 U x rdr + / pvu' U r rdr 

Jo Jo 

/•OO fOO 

+ / pv'v'V r rdr+ / pV w'w'dr 
Jo Jo 

1 r°° — „ 

“ / {U r ) 2 rdr 


Re 


(30) 


’o 


1 

3Re 


M'*(V • v') x + l/*(V • v') r 


Similarly, the energy equation for the large-scale struc- 
ture component reduces to; 


w'* 

+ — (V- «')* + «; 

where Q = (|w/| 2 + |n / | 2 -|-|u/| 2 )/2, cc denotes the complex 
conjugate of the preceding complex terms, and ht denotes 
higher order terms containing the product of three or more 
disturbance variables. 

When the shape assumption given by equation (8) is ap- 
plied in equation (28), the axial derivatives will decompose 


d r°° 

dxl (WQMr = 




J ^'*£4 + v'*p' r H — + ccj rdr 

J (u'*pu U x + 1 i'*pv U r + ccj rdr 

_ pw'*w ' — \ 

pv'*v' V r + V + cc I rdr 
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1 

Re 


•y w 

+-o-K + M,) + -t(»' - 2 <) + cc 


(31) 


rdr 


and 


( 34 ) 


+ 



• v') x + t/*(V • v') r 


w 


(V-v 1 ) 


cc 


rdr 


All the integral terms in equations (32) to (34) are given 
in the Appendix. Note that unlike in the previous paper, 4 
all the energy transfer integrals appearing in equation (32), 
represented by I m wt, appear identically in equation (33) 
with opposite sign, as expected. 


In this equation, terms to second order are retained and 
those terms involving the second derivative with respect 
to x are neglected. For terms higher than second order, 
the time-average is zero when using the shape assumption, 
equation (8). 

The physical interpretation of the terms appearing in the 
energy equations is obvious. For the mean flow equation, 
the first term on the left side is the mean flow advection 
of the mean flow kinetic energy, the second term involves 
energy transfer. The first four terms on the right-hand side 
of (30) govern the energy transfer from the mean flow to 
the coherent structure, and the last term is the viscous dis- 
sipation of the mean flow energy. As for the large-scale 
equation (31); the left-hand side is the mean flow advec- 
tion of the large-scale kinetic energy. The first term on the 
right-hand side is the work done by the large-scale pressure 
gradients, the next two terms are the energy transfer from 
the mean flow to the large-scale structures, and the last two 
terms are the large-scale energy dissipation. 

When substituting the mode definitions into equation 
(30) for the mean energy, we will assume that the distur- 
bances are real. Hence, each disturbance variable is rep- 
resented by equation (8) plus its complex conjugate. The 
time-averaging process results in the following ordinary 
differential equation. 

dddl^ + d = a2 1 

dx d9 dx L J Re 

In this analysis, the mean flow is characterized by the mo- 
mentum thickness rather than by the axial distance. Hence, 
the integrals in the energy equations are dependent on 6 
instead of x. 

For both the large-scale energy equation (31) and the 
phase equation, the disturbance variables are used in their 
complex form as given in equation (8). The resulting ordi- 
nary differential equations are: 


d 

dx 


[■ A 2 I a 


dA ? J _ A 2 dl toT 

dx Zll wpla 


(33) 


= A 2 



Ipt 


1 

Re 


d'/ntl 



NUMERICAL RESULTS 

The calculations were performed for three, perfectly ex- 
panded jets with exit Mach numbers of 1.5, 1.8, and 2.1. 
In all cases, the jets were cold with the jet exit total tem- 
perature equal to the ambient temperature. The mean flow 
was computed using the parabolized, boundary-layer pro- 
cedure of Dahl & Morris 13 with the addition that the radial 
velocity profile was computed from the mean continuity 
equation. This procedure generated mean flow profiles that 
smoothly transition from the initial region to the fully de- 
veloped region allowing complete parameterization by the 
momentum thickness, 6, as required by the nonlinear anal- 
ysis. 

The mean flow was used in the Linearized Euler Equa- 
tion (LEE) solver described in a previous section. The 
time-domain solver was run until initial transients ended. 
The unsteady disturbance flow field was then stored at con- 
stant time intervals creating time histories of all the com- 
puted unsteady flow variables at each grid point of the flow 
field. The time interval corresponds to 1 / 8th of the period 
for the frequency used for the inflow disturbance. In our 
calculations here, the frequency corresponds to a Strouhal 
number fD/Uj = 0.2 and n = 1 for the azimuthal mode 
number. An example of the computed pressure disturbance 
field is shown in Figure 1. The pressure disturbance con- 
tours exhibit a growth and decay behavior within the jet and 
a radiation directivity characteristic outside of the jet. The 
pressure disturbance also decays with distance away from 
the jet. These results are typical for all three computed jets. 

The set of disturbance flow fields from the LEE cal- 
culations were post-processed to provide the data for the 
nonlinear integral calculations. After taking the Fourier 
Transform in time to obtain the complex amplitude of each 
unsteady flow variable at each grid point, the data were 
normalized according to equation (9) to obtain the radial 
shape functions for each variable. The radial shape func- 
tions are then used to compute the integrals that appear in 
equations (32) to (34), the set of coupled ordinary differ- 
ential equations for the mean flow momentum thickness 9, 
the disturbance variable amplitude A 2 , and the disturbance 
variable axial phase function ip. Since the integrals are con- 
sidered to be functions of 9 , interpolations were made of 
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these integral functions during the solution of the differen- 
tial equations. The differential equations were manipulated 
to form three first-order, coupled, initial-value differential 
equations that were solved by implementing a degree- 1 
Taylor series scheme with fixed point iteration. 14 Only the 
choice of the initial amplitude was needed to start the inte- 
gration at x = 0. The initial value for 9 was given by the 
mean flow calculations and the initial value for the phase 
function ip was set to zero. Any nonzero initial value for 
ip merely offset all the results for ip by that initial nonzero 
value. 

Effect of Re 

The differential equations contain a Reynolds num- 
ber that comes from the normalized Navier-Stokes equa- 
tions. Obviously the Euler equations are inviscid, hence 
the Reynolds number must be provided by the underlying 
mean flow. Solutions for the amplitude function are shown 
for the Mj = 1.8 jet in Figure 2 for two different ini- 
tial amplitudes, A% = 1CT 7 and A ° = 10“ 3 . Given that 
we consider the jets to be high-Reynolds number turbu- 
lent jets, the high constant Reynolds number results show 
growth but very little decay for the amplitude. If we con- 
sider the viscosity to be a turbulent viscosity, as in Tam & 
Chen, 15 and choose a smaller constant Reynolds number, 
we see that significant ’early’ decay does not begin until 
the Reynolds number is getting relatively small. To avoid 
the arbitrary choice of a Reynolds number, we follow our 
previous procedure by choosing a local variable Reynolds 
number obtained from the mean flow calculations. 4 This 
is a simple way to mimic the increased effects of fine- 
scale turbulence on dissipating the energy in the large-scale 
turbulence as the large structure moves downstream. The 
amplitude results using the variable Re approach is also 
shown in Figure 2. 

Local energy integrals 

The four integral terms in equation (32) govern the ad- 
vection, transfer, and dissipation of the mean kinetic en- 
ergy. Using the data computed by the LEE code, these 


Fig. 2 Amplitude function comparisons showing effects of 
Reynolds number. Mj = 1.8; fD/Uj = 0.2; n = 1. 




Fig. 3 Energy integrals from the mean flow energy equation 
(32). fD/Uj = 0.2; n = 1. 


integrals are computed using the trapezoidal rule and the re- 
sults are plotted in Figure 3 for the three jet Mach numbers. 
Due to the transient waves present near the inflow bound- 
ary, the integral calculations were started at 9/D = 0.035. 
As functions of 0 , the two integrals that depend solely on 
the mean flow, I am and I m d , are nearly identical for the 
three jets. The integral I am represents the advection of the 
mean flow kinetic energy. The rate of change of I am is 
initially constant and then gradually decreases. The energy 
moves downstream as the jet spreads and then the advection 
slows downstream of the potential core. The I rn( [ integral 
represents viscous dissipation of the mean kinetic energy. 
It starts out initially high due to the high radial gradients of 
the axial mean velocity. It steadily decreases in value as the 
flow proceeds downstream. 

The I mw t integral governs the main transfer of energy 
from the mean flow to the large-scale structures. It is en- 
tirely negative as the jet evolves, indicating that energy 
is only transferred from the mean flow to the large-scale 
structures and not the reverse. Its magnitude is largest 
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e/D 


Fig. 4 Comparison of energy transfer integrals. I mw i is 
dominant integral. I mw t is Imwi plus second order integrals. 

for the lowest Mach number and it decreases as the Mach 
number increases becoming a less efficient energy trans- 
fer process. In the current derivation, I mw t consists of the 
sum of six terms, the largest, I mw 1 , has the classical form 
for the velocity-shear driven generation of the large-scale 
structures (See Appendix). To various degrees, the addi- 
tional terms add effects of compressibility and non-parallel 
flow to the energy transfer process. As can be seen in 
Figure 4, these additional integrals add little to the total 
energy transfer integral in these cold, slowly diverging jets. 
The waviness seen in the initial part of the plots of I mw t 
in Figure 3 are remnants of the problems with the inflow 
boundary condition. 

The only way to possibly transfer energy from the large- 
scale structures back to the mean flow is via the I at integral. 
The general behavior of I at is to increase to a maximum 
value and then decrease downstream. The initial positive 
gradient is indicative of energy transfer towards the large- 
scale structures. After the peak and the gradient becomes 
negative, then some energy can be transferred from the 
large-scale structures to the mean flow. Both the magni- 
tude and the gradients each side of the peak increase as the 
Mach number increased. 

Figure 5 shows four of the dominant integral terms in 
equation (33) governing the kinetic energy in the large- 
scale structures. The integral for the advection of the large- 
scale kinetic energy, I aw , increases as the energy flows in 
from the mean flow. Once dissipating mechanisms take 
over, the integral decreases in value. The lower values of 
I aw as Mach number increases corresponds to the lower 
levels of incoming energy from the higher mean flows. 
The large-scale energy decreases as I w d , the dissipation 
integral, increases. The remaining two integrals shown 
represent work done by the large-scale structure pressure 
gradients; I p i represents work done by the axial pressure 
gradients and I pt represents work done by both the radial 
and the azimuthal pressure gradients. Both integrals are 



e/D e/D 


Fig. 5 Energy integrals from the large-scale structure energy 
equation (33). fD/Uj = 0.2; n = 1. 

relatively small in magnitude. These integrals are the only 
ones using the pressure disturbance and their behavior at 
the computed flow field exit indicates the presence of re- 
flections from the outflow boundary. 

Nonlinear amplitude 

The results for the computed amplitude function for the 
jet with Mach number 1.8 are shown in Figure 6. The ini- 
tial amplitude A 2 varies from 10~ 7 to 10 -1 . Part (a) shows 
the absolute value of the amplitude and part (b) shows the 
amplitude normalized by the initial amplitude. Notice that 
when the amplitude is normalized, the A 2 = 10 1 and 
1 0” 6 plots are basically identical. This indicates a linear 
response. The amplitude function changed by the same 
amount as the change in the initial amplitude. The nor- 
malized plot for A 2 = 10~ 5 begins to deviate from the 
previous normalized plots at about x/D = 10. In part (a), 
it is seen that the A 2 = 10 -5 plot is exceeding approxi- 
mately 10 -3 at that location. This pattern holds for higher 
initial amplitudes indicating an amplitude threshold above 
which nonlinear effects take place. At the initial amplitude 
of 10 -1 , growth has nearly ceased and single frequency 
saturation starts to occur. 16 

Using initial amplitudes of 10 2 and 10 to generate 
nonlinear and linear responses, respectively, for the ampli- 
tude function. Figure 7 shows comparisons for the three jet 
Mach numbers. The amplitude of the linear response grows 
in amplitude as the Mach number increases. As the ini- 
tial amplitude increases and nonlinear effects become more 
important, the Mach number has very little effect on the 
maximum large-scale structure amplitude. 

The energy integral equations also compute the change 
in mean flow momentum thickness and large-scale struc- 
ture phase function. Results are shown in Figure 8. Part (a) 
shows that the momentum thickness remains unchanged for 
small initial amplitudes. It is through the nonlinear process 
that the momentum thickness increases at the higher am- 
plitudes until saturation. In part (b), the momentum thick- 
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x/D 


Fig. 6 Amplitude function comparisons with different initial 
amplitudes, (a) Absolute amplitude, (b) Amplitude normal- 
ized by initial amplitude. Mj = 1.8; fD/Uj = 0.2; n = 1. 



x/D 


Fig. 7 Amplitude function comparisons between nonlinear 
and linear calculations at different Mach numbers. Non- 
linear: A 2 = 10 -2 . Linear: A 2 = 10 7 . fD/Uj = 0.2; 
n = 1. 

ness decreases downstream as the Mach number increases 
where initially the behavior of the momentum thickness 
is similar for all three jets. Finally, the phase function is 
plotted in part (c) as a phase speed for the large-scale struc- 
ture. Nonlinearity has the effect of a rapid increase in phase 
speed near the nozzle exit compared to the gradual change 
in phase speed seen in the linear case. 

CONCLUDING REMARKS 

The Linearized Euler Equation method, with given non- 
parallel mean flow, has been used to compute the large- 
scale turbulent structures. Normalized radial functions 
were derived from this data to be used in the energy integral 
method that describes the weakly nonlinear interactions of 
a growing large-scale structure with the mean flow. Below 
some initial amplitude threshold, the evolving disturbances 
behave linearly. Above that threshold, the amplitude of the 
disturbance peaks and decays earlier in axial location com- 
pared to the linear response. These results were also shown 
to be highly dependent on the dissipation model used for 
the large-scale structures. 
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